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Abstract 

Fermi acceleration is the process of energy transfer from massive objects in 
slow motion to light objects that move fast. The model for such process is a 
time-dependent Hamiltonian system. As the parameters of the system change 
with time, the energy is no longer conserved, which makes the acceleration 
possible. One of the main problems is how to generate a sustained and robust 
energy growth. We show that the non-ergodicity of any chaotic Hamiltonian 
system must universally lead to the exponential growth of energy at a slow 
periodic variation of parameters. We build a model for this process in terms of 
a Geometric Brownian Motion with a positive drift and relate it to the entropy 
increase. 
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1 Introduction 


Consider a family of Hamiltonians H{p, g; r), where r denotes a set of param¬ 
eters changing slowly with time. The equations of motion are 


q 


dH 

dp 


{p,q;T), 



{p,q;r). 


( 1 ) 


As parameters change, the energy E 
system: 

E = 


= H{p, q; t) is not preserved by the 


dH . 


( 2 ) 
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Our goal in this paper is twofold; we show that under very general conditions 
periodic oscillations of parameters lead to the exponential energy increase, and 
we investigate the main features of this process. 

Physically, model Q accounts for the energy transfer from massive objects 
in slow motion to light objects performing a fast motion. The (p, q) variables 
are momenta and coordinates that correspond to the fast degrees of freedom, 
while the parameters r describe the effect of the slow degrees of freedom (we 
assume that they correspond to a certain massive object, i.e., their evolution 
is not influenced by the fast dynamics). It is known since the Fermi’s work 
on cosmic rays pQ that acceleration is possible in this setting, see also [2]. In 
addition to cosmic rays studies [3] , the interest to the problem is also motivated 
by plasma confinement [1] and nuclear fission [5]. 

Previous research has been most often focused on time-dependent billiard 
models. Billiard is a dynamical system corresponding to a particle which 
moves inertially inside a bounded domain and, upon hitting the boundary 
of the domain, reflects elastically. In time-dependent billiards, the boundary 
moves according to a certain given law which is not affected by the collisions 
with the particle (i.e., the billiard wall is infinitely heavy). The wall motion 
changes the particle reflection law, so the particle’s kinetic energy is no longer 
conserved at the collisions. Still, this does not lead to an acceleration in 
the one-dimensional case (where the particle moves inside an interval with 
slowly and periodically moving end points). Here, the kinetic energy stays 
bounded for all times m 0 0 m, provided the law of the end points’ motion 
is smooth enough. As was discovered in [lOI lllj . in the two-dimensional case 
the situation changes. When the domain’s boundary slowly moves so that the 
dynamics inside the corresponding frozen billiard is chaotic the time-dependent 
billiard produces accelerating particles, as can be seen from various numerical 
experiments [la [HI la [la [Ml Ea EE]. The theory of this phenomenon was 
even earlier proposed by Jarzynski who derived a universal model describing 
the evolution of the distribution of energies in a strongly chaotic billiard with 
moving walls m- In these examples, the observed acceleration was relatively 
slow, giving at most polynomial in time growth of the energy averaged over a 
uniform ensemble of initial condition. Such regime can be easily destroyed by 
a small dissipation I20l|2l|. 

It was shown in |221 [23] that the main limitation for fast energy growth is 
imposed by the existence of an Anosov-Kasuga adiabatic invariant |24l [25l [26l 
EH EH]. This adiabatic invariant is a function of the particle energy and the 
slow time r. It follows from the Anosov averaging theorem 1211E9] that the 
change of its logarithm over the period of the billiard wall oscillations tends to 
zero for the majority of initial conditions as the particle energy grows. This 
means that the rate of possible acceleration must slow down with the energy 
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increase. The logarithm of Anosov-Kasuga invariant coincides with the Hertz 
entropy of the frozen billiard [30]. Therefore, its approximate preservation is 
in full agreement with the basic physical intuition: the entropy of an ergodic 
system does not grow at slow (i.e., adiabatic) changes of system’s parameters. 

Recently, various examples of time-dependent billiards with exponential 
energy growth were produced |3ll \22l EZ] [33l [3l|. In these works, the billiard 
wall moves in such a way that the frozen billiard loses its ergodicity (i.e., it 
has more than one ergodic component) at least at the part of the period of 
the boundary oscillations. These works ergodicity violation, but producing 
the same effect: The Anosov averaging theorem does not hold for non-ergodic 
systems, so the Anosov-Kasuga invariant no longer exists, and the energy grows 
unimpeded with each period of the billiard wall oscillations. 

In this paper we depart from the billiard setting, and investigate a general 
question: How can an adiabatic periodic variation of parameters lead to a sus¬ 
tained energy growth? The Anosov-Kasuga adiabatic invariant is not billiard 
specihc, so it imposes restrictions to the energy growth in any Hamiltonian 
system with slowly changing parameters if the frozen dynamics is ergodic on 
every energy level. However, apart from special classes of systems, such as 
geodesic flows and billiards, a typical Hamiltonian system is not ergodic. 

We demonstrate that the non-ergodicity of a chaotic Hamiltonian system 
must universally lead to the exponential growth of energy at a slow periodic 
oscillation of parameters. The key mechanism is the following: A non-ergodic 
Hamiltonian system has regions of chaotic dynamics in the phase space, which 
coexist with stability islands where dynamics is nearly integrable (quasiperi- 
odic). Adiabatic changes of parameters lead to transitions between these re¬ 
gions. Different initial conditions give rise to different itineraries of these tran¬ 
sitions, and different itineraries give different values of the energy gain/loss 
per period of the parameters oscillation. We introduce an analogue of Hertz 
entropy for non-ergodic systems and show that on average over all possible 
itineraries the entropy linearly increases after each period. This yields the 
exponential energy growth. 

1.1 Assumptions and claims 

Consider a family of Hamiltonians H{p, q-, r) where q,p ^ and r is a set of 
parameters. We assume that 

Al the Hamiltonians in the family are homogeneous at each frozen r, 

A2 the parameters r change periodically with time, 

A3 the Hamiltonian H(p,q]T) has more than one ergodic component (on 
each energy level) at each frozen value of r for at least part of the period. 
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Our hypothesis Al means that each Hamiltonian is invariant with respect 
to energy scaling, so the dynamics on each energy level is the same. A typical 
example is the motion in a homogeneous polynomial potential, see e.g. (22). 
Another example is the Boltzmann gas of hard spheres. Note that assumption 
Al is not restrictive. We may start with an arbitrary system, but since we 
study the process of an unbounded energy growth, we must consider the system 
in the high-energy limit, and it is most natural to assume that the system 
becomes invariant with respect to the energy scaling in this limit. For example, 
if we start with a general polynomial potential, only the highest order terms 
will be most relevant at high energies, i.e., the potential will effectively become 
homogeneous. For the motion in a potential confined to a bounded domain 
D (i.e., the potential which is finite inside D and infinite on the boundary of 
D) the high-energy limit is the billiard in D [35], which is also a homogeneous 
system. In essence, our homogeneity assumption means that the external 
forcing is applied at the highest order terms, i.e., the forcing is of the same 
order as the energy: 

Unless we consider relativistic particles, the values of q and p grow with 
the energy while r remains bounded. Thus, at large energies, the variables 
(p, q) are fast and the parameters r are varying slowly. Another slow variable 
is the energy E = H{p(t), q{t);T{t)): since dH/dr ~ H, the speed of change 
of InF^ is comparable with f. 

A proper way to understand the evolution of slow variables is to somehow 
average the system over the fast variables. Such model reduction goes naturally 
when the frozen system is ergodic with respect to the Liouville measure, i.e., 
Lebesgue measure restricted to the constant energy level; 


pLL{E,T;dpdq) = 6{E - H{p,q;T))dpdq. 

In this case Anosov averaging theorem |24l [29] guarantees that the evolution 
of the energy can be described by averaging equation ([^ over that is 

I ^{p,Q-,T)pLiE,T;dpdq) 

j PL{E,T\dpdq) 

with good accuracy and for a large set of initial conditions, see details in 
Section In the non-ergodic case there is no unique measure the averaging 
over which could provide a description of the slow evolution of energy which 
would be valid for the majority of initial conditions simultaneously. Indeed, 
below we propose a model where the energy change is governed by equations 
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obtained by averaging over certain randomly chosen ergodic components, exact 
choice of which depends on the initial conditions. 

The main approximation principle we use in our approach is that we as¬ 
sume that at every value of the slow variables the motion of the fast variables 
immediately ergodises. One can think of this as at each value of r and E the 
phase point is uniformly distributed over a certain ergodic component of the 
fast system. We stress that we view this only as an approximation, i.e., our 
assumption is that the time averages over intervals long in terms of the fast 
variables and short in terms of the slow variables coincide with averages over a 
certain ergodic component of the fast system with a sufficiently good accuracy 
for a sufficiently large set of initial conditions. 

In this approximation, at each value of r and E the evolution of E is defined 


by 




(4) 


where the averaging is done over a certain ergodic invariant measure of the 
frozen system on the corresponding energy level. Moreover, the fast ergodis- 
ation assumption also implies that the probability to switch in this equation 
from one ergodic component to another (as r changes) depends only on the 
pair of ergodic components involved. Thus, we have a time-dependent (peri¬ 
odic) Markov process on the set of ergodic components of the frozen system. 
Without losing much precision one can discretise time and also approximate 
the set of ergodic invariant measures by a finite set, so one obtains a process 
of hopping over a finite set of states. Different initial conditions correspond to 
different realisations of this process, i.e., to different time-dependent families of 
ergodic invariant measures; we call such a family an averaging protocol. Each 
protocol produces its own averaged equation for the evolution of energy. 

Similar construction was described in [3H |36]. The validity of such ap¬ 
proach was numerically checked in [221 [551155] for various examples of time- 
dependent billiards. In these examples, the structure of ergodic components 
was known, so the transition probabilities for the hopping process were com¬ 
putable. In a general situation (e.g. in the example we consider in Section 
there is no hope of knowing the structure of the set of ergodic components, 
nor computing the transition probabilities. Therefore, we formulate the main 
properties of the resulting model, which are independent of its particular de¬ 
tails. 


Multiplicative law for the energy evolution. The first conclusion we 
make is that after each period T of the parameters’ oscillations the energy is 
multiplied to a random factor: 

En+l — E^Qni (5) 
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where = H[p,mT{nT)) is the energy after n periods T. The multiplicative 
character of law (IM is due to our assumption ([^ which implies that E ^ E 
in 0 no matter which ergodic component is chosen to average over. The 
energy gain factor depends on the particular averaging protocol. As the 
sequence Cn is formed via a Markov hopping process which typically exhibits an 
exponential decay of correlations, the correlations between consecutive values 
of Cn niust rapidly decay. For simplicity we just assume that the factors Cn 
are independent random variables; they are also identically distributed (as the 
system we consider depends on time periodically). 

Exponential growth of energy. Model ([^ describes a random walk for 
In En- By the Law of Large Numbers, for a typical realisation of the random 
walk (i.e., for a typical initial condition) we must have 

lim — InEn = Eln(). 

n^+oo n 

In Section we introduce a notion of entropy for our system and show that 
the preservation of volumes in the {p, q) space by system Q implies that the 
random process ([^ cannot decrease the entropy. In the ergodic case the en¬ 
tropy would be preserved, but in the general case of an adiabatically perturbed 
homogeneous Hamiltonian system with several distinct ergodic components in 
each energy level we expect that the entropy increases to a non-zero quantity 
with each period T. 

This is equivalent to the claim that the multiplicative random walk model 
([^ has a positive bias: 

p = ElnC > 0. (6) 

Hence, the energy grows exponentially both for typical initial conditions and on 
average. As we mentioned, the non-ergodicity plays an important role here: In 
the ergodic case the bias p vanishes and model ([^ becomes invalid (the energy 
grows at most polynomially [221123| ). 

Note that since the factors Cn in @ are assumed to be independent, the 
energy averaged over all initial conditions grows by the law 

EEn+l = (ECn) EEn. 

This corresponds to the exponential growth with the rate 

— lnEE„ = = InE^. 

n 

As p'^ > p, the averaged energy growth is faster than the energy growth for 
a typical conditions, which means that a small minority of realisations far 
outperform the rest. 
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Log-Gaussian distribution. The positivity of the rate p is the universal 
characteristics of the energy growth process in adiabatically perturbed non- 
ergodic systems. The distribution of energies formed at this process is not 
universal, as billiard examples of [22l [33l show. As In En is, according 
to ([^, the sum of independent bounded random variables, the logarithm of 
the energy becomes, after a proper scaling, normally distributed at large n. 
This does not mean that the energy itself acquires a Log-Gaussian distribution 
(since the scaling and taking a logarithm do not commute). However, in the 
example we consider in this paper the energy distribution is, in fact, close to 
Log-Gaussian, right after the first period (see Figs. 7p ). 

We believe that this must be typical for the case where the frozen Hamilto¬ 
nian system has many ergodic components (at least for a part of the period T), 
so a typical trajectory will hop many times between the components during 
the period, see Figure for an illustration. By Q and (§ , if for a certain 
time interval the change of energy along the orbit is determined by a given 
ergodic component, the energy at the end of this interval equals to the energy 
at the beginning of the interval times a certain factor; the product of these 
factors over the period gives the multiplier from Q. Thus, if the number 
of jumps between the components is large, then occurs to be a product of a 
large number of random factors, each of which must be close to 1 (as C must 
remain bounded). This means that the logarithmic gain In ^ at the end of the 
period is in the Gentral Limit Theorem regime, i.e., the distribution of ln(( is 
close to a Gaussian. Since the sum of independent Gaussian variables is also 
Gaussian, we conclude that the distribution of the scaled energy after 

n periods is close to Log-Gaussian; moreover, we can estimate parameters of 
this distribution: 

^([nEn) ^'y^{np,y/na) (7) 


where ^(a, 6) denotes the normal distribution with the mean a and standard 
deviation b, and 

0-2 = E(lnC - p)^ = 2(/9+- p). (8) 


We note that except for specially prepared examples (like above mentioned 
billiards) it is impossible to theoretically compute the rates p and p'^; however, 
we believe that relation Q between the energy growth rates and the width 
of the distribution of the logarithmic gain must hold for a large class of (not 
specially prepared) systems. 

This particular class of multiplicative random processes, where the log¬ 
arithmic gain is Gaussian, is called geometric Brownian motion (GBM). It 
emerges in various applications, notably in finance, and is well studied, see 
e.g. [371 EE] and references therein. 
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2 Entropy growth 

Now we give a more detailed description of the model, introduce the notion of 
entropy, and demonstrate the validity of Q. We call a Hamiltonian H{p,q) 
homogeneous if for any E > {) there exists a coordinate transformation 
that keeps the system the same, sends the energy level H = 1 to H = E, and 
has a constant Jacobian 

J{E) = E^, a> 0. 

We assume that the positive energy levels are compact, so 

V{E) = V{l)J{E) = (9) 

where V(E) is the volume of the {p, g')-space between the energy levels H = E 
and H = 0. Thus, we can label the points in the phase space (p, q) by the 
coordinates {x,E) where E = H{p,q) is the energy and x = g) is the 

projection to the energy level = 1, see Fig. 



Figure 1: Illustration to the (x, ii^)-coordinates. 


We consider a r-dependent family of homogeneous Hamiltonians H(p, q; r), 

so the transformation $ and the volume V are now also functions of r. We 
d (dH\ 

assume that — InJJ = ( j /H does not depend on energy.This means that 


dr 

^ scales in the same way as H and condition 


dH 
dr 

the notation 


holds. We also introduce 


G{x,t) := -^\nH{q,p;T). (10) 

Let us now allow the parameters to change periodically with time. We 
assume that at large energies the motion in the frozen system (that corresponds 
to a fixed value of r) is fast, i.e., the variables x change, at large energies, much 
faster than r does. By (|^, the speed of change of InF is comparable with f. 







Thus, we have a slow-fast system, with fast variables x and slow variables r 
and Ini?. 

If the frozen system is ergodic on any energy level with respect to the Liou- 
ville measure /r^, we can apply Anosov theorem [2^. This theorem guarantees, 
that for any given finite number of periods and any given accnracy, if the initial 
value of energy is taken large enongh, then the evolntion of the logarithm of 
energy is, within this accuracy, described on the chosen time interval by the 
averaged equation 


^InE = 
dt 


dH 

-^{p,q,r)6{E - H{p,q,T))dpdq 
E [ 5{E — H{p,q,T))dpdq 


( 11 ) 


for the majority of initial conditions (i.e., for all initial conditions except for a 
set of small measure). By changing variables {p,q) to {E,x = ^f}{p,q)), we 
reduce this equation to 


d . ^ 1 

— In i? = -p— 
dt J dx 


G{x, T)dx 


( 12 ) 


where G is defined by (10); the variable x rnns the level H = 1. 

We note that no mixing is required, i.e., one should not think of a fast 
relaxation to the Liouville measure at every moment of the slow time. Just 
the measure of the set of initial conditions for which the evolution of energy 
deviates noticeably from that given by this averaged eqnation is small (the 
measure of the bad set of initial conditions can be made as small as we want 
by taking the initial energy value large enough). 



Figure 2: Illustration of the jumping process between macro-states. Initial points 
will have different itineraries visiting distinct ergodic components during the period. 


Note that equation 0 implies 


fL 

dt 


V{E,t) 


dV ■ dV . ^ 
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Indeed, the volume under the energy level H = E is given by 


V{E,t)= f dpdq = f 9{E — H{p,q,T))dpdq 

JH(p,q,T)<E J 

where 9 is the Heaviside function. Since 9' = <5, it follows that 

/ -^iP,q,T)6{E - H{p,q,T))dpdq, 


dV 

Ih 


= j 6{E - H{p,q,T))dpdq, 


so 


^y(£.r)=0 

by ([IT]) (cf. [25l |26l [271 ESI [3^ ). This means that energy changes periodically 
with r (to keep V{E,t) constant). As the behaviour of InE in the original 
system is close to that given by (11), we conclude, using ([^, that for the large 
majority of initial conditions In V{E, r) stays close to its initial value for the 
given in advance number of periods, Inii^ stays bounded and, at the end of 
each period, returns close to its initial value InSoj he-, E/Eq returns close to 1 
at the end of each period. This shows that if the frozen system is ergodic, then 
there can be no significant energy growth (nor decrease) for a long interval of 
time and a large set of initial conditions. 


In the non-erdodic case this conclusion is no longer valid. Moreover, it is 
no longer true that the evolution of energy is approximately the same for the 
majority of initial conditions. Maximum we can hope, in analogy with the 
Anosov averaging theorem, that the slow evolution of the energy is described 
(for a given number of slow periods, with a prescribed accuracy, and for a large 
set of initial condition) by the following equation (cf. (dl)): 

^ln£'= G{x,T)pr{dx'^ f (13) 

where pr is, at each value of r, a certain ergodic measure (invariant with 
respect to the frozen system) on the energy level H = 1. Moreover, as we are 
looking only for a finite accuracy approximation, we can replace the measures 
Pr by one of a finite set of measures 

dx 

/i?=i 


Pit = ii{T,x) 
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where at every given r, are characteristic functions of certain open sets 
which form a partition of the x-space; we also may assume that the dependence 
of on T is piece-wise constant. 

We call the measures m (or their support sets) macro-states. The orbit of 
the system hops, as r changes, between different macro-states. As we men¬ 
tioned in the previous Section, a part of our fast ergodisation principle is that 
the corresponding sequences of macro-states (the averaging protocols) are re¬ 
alisations of a r-dependent Markov process. For a typical Markov process, for 
any given distribution of initial conditions over the macro-states, after a few pe¬ 
riods of oscillation of r the distribution will converge to a certain time-periodic 
distribution. Thus, independently of the given initial density of phase points 
in the (x, i?)-space, we may assume that after a few periods T the dynamics 
will form the same piecewise constant density in the x-space at the beginning 
and the end of each period; this equilibrium density will take constant values 
on macro-states. 

The next question is to investigate the dynamics of the distribution of en¬ 
ergies. Note that in our construction we have a finite set of averaging protocols 
over one period T of oscillations of r (the number K of the protocols depends 
on the accuracy of the approximation we want to achieve, but it does not 
depends on the energy, i.e., our model provides the same accuracy of approxi¬ 
mation on a given number of periods for all sufficiently large values of energy). 
Let us split the space of initial conditions x into meso-cells Mi ,..., Mk that 
give rise to the distinct averaging protocols. By construction, each cell is a 
subset of one of the macro-states that exist at the beginning of the period, 
and the image of the cell by the flow of the system after the period T is also a 
subset of one of these macro-states. For each cell the majority of initial condi¬ 
tions gives rise to the same energy evolution (given by equation p^) over the 
period T, while for initial conditions from different cells the values of energy 
gain or loss will be different. 

Let Eq and Ei be two sufficiently large values of energy. Choose a given 


cell Mfc (fc = 1,... ,K). By equation (13), if the orbits with initial conditions 
E = Eq, X ^ Mj^ move to the energy level E = Eq = e^’^Eo after the period 
r, then likewise, the orbits with initial conditions E = Ei, x € move to 
the level E = Ei = e^'^Ei. Now note that the non-averaged system preserves 
volume in the {p, g)-space. Therefore, it follows that the volume occupied by 
the points with x G between the levels E = Eq and E = Ei equals to 
the volume occupied by the points with x G between the levels E = Eq 
and E = El, where denotes the image of the cell by the flow of the 
(non-averaged) system after the period T. This gives 


viMk)iViEi) - ViEo)) = v{Mk)iV{Ei) - V{Eq)) 
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where v is the volume in the x-space (the level H = 1), and V(E) is the volume 
under the energy level H = E (all volumes here are computed at the beginning 
of the period T, or at the end, which is the same). By plugging © into this 
formula, we obtain the following relation between the energy growth rate Xk at 
the cell Mfc and the volumes of the cell and its image Mf^ after one period: 


aXk = In 


v{Mk) 

viMkY 


(14) 


We also normalise v so that the total volume of the x-space at the beginning 
of the period equals to 1, that is. 


^u(Mfc) = ^u(Mfc) = l. 


(15) 


Given a distribution v of the points in the (x, E)-space, we define the 
corresponding entropy as 


S= In 


V{E,t) 
V{1 


^^v{dE,dx) = a [ In E n(dE,dx), 

,V J 


(16) 


see Q. In our model, the energy undergoes the same evolution for the majority 
of the orbits that start at the cell — for all of them In E gets an increment 
Xk- Therefore, the change of the entropy over the period equals to 


AS = a'^ XkVk 


where Uk is the probability to be at the cell Mk at the beginning of the period 

^(Mfc) \ 


By (14), we obtain 


^ In 

k 


YMk)J 


Vk- 


(17) 


The result depends on the distribution v. However, as we mentioned, in 
our model any initial density in the x-space converges to an equilibrium one, 
constant on macro-states. For example, there can be a situation where at some 
value of r during the period the frozen system appears to be completely chaotic, 
i.e., at this value of r there is only one macro-state, hence, the corresponding 
equilibrium density is uniform. In this case, we choose this value of time to be 
the beginning of the period, so Vk = v{Mk) in © (recall that v is the scaled 
volume so that (15) holds). Note that we do not insist on ergodicity or mixing 
at this value of r. All our formulas can be true only within certain, chosen 
in advance finite accuracy, so we just need that the uniform density would be 
close to our equilibrium density with the given accuracy. This seems to be 
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the situation we have in the example considered in Section]^ the numerically 
obtained phase portrait at the beginning of the period seems fairly uniform. 
Thus, if we have only one macro-state at the beginning of the period, then 


AS 


'^Vk In 

k 


' 1 ^ 

Vk 


where Vk = v{Mk), Vk = v{Mk). As Ylvk = 1, and the geometric mean is 
smaller than the arithmetic mean, we have 



= 0 


(see (EH)- Note that AS* can be zero only if Vk = Vk for all k. By ( |14[ ), this 
would mean no energy change for each averaging protocol, which seems to be 
a quite degenerate phenomenon for the non-ergodic case. Therefore, unless we 
are in some degenerate situation, the entropy must acquire the same positive 
increment at each period. By the definition of entropy (16), the linear growth 
of entropy means an exponential energy growth. 


The same conclusion holds in the general case where we have more than 
one macro-state at the beginning of the period. In this case we have a certain 
equilibrium density in the x-space, which produces the probabilities Vk of being 
at the cell Mk, the same at the beginning of each period. Recall that the 
equilibrium density is constant on a macro-state, and that each cell Mk is a 
subset of a certain macro-state, so the equilibrium density is uniform on Mk- 
Therefore, the equilibrium density is given by 

u{x) = Vk/v{Mk) for X G Mk, (18) 

where v is the volume in the x-space. The points that start at a cell Mj come, 
after time T, to the image cell Mj. It is a part of a certain macro-state, and 
our fast ergodisation principle (see Section tells us that the probability of 
a point from Mj to be anywhere in this macro-state is the same within the 
accuracy of our approximation. Therefore, the probability Qkj to start at the 
cell Mj and, after time T, land at the cell Mk, equals to v{Mj n Mk)/v{Mj). 
If we denote Vk = v{Mk), Vj = v{Mj), we obtain 

Vk — ^ ^ QkjVj ( 19 ) 

j 

(as v{Mk) = Ylj v{Mj n Mk)). Recall that the probabilities i^k are the same at 
the beginning of each period, so we also have 

k'k ^ Qkj k'j . ( 20 ) 

j 
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By the change of entropy over the period S is given by 


AS" = Uk In 

k 


Vk 


= s-s, 


( 21 ) 


where 


S = ^ i/fc In 


k'k 


— / i'{x)\n.v{x)dx and £ = Vj In 

j 




i.e., £ is the Shannon-Gibbs entropy of the equilibrium density (18), and £ is 
the Shannon-Gibbs entropy of the density ^'(x) = ^ at x G Mj. The density u 
is the image of by the slow evolution over the oscillation period, just before 
it is replaced by u due to the fast ergodisation. The ergodisation must increase 
the entropy, so AS > 0. 

Indeed, let jkj = QkjVj/i'k- Since Yhj Ikj = 1 (see (20)), and the geometric 
mean is less than arithmetic mean, we have 


-ikj Hvj/T^j )=In (^) j 



which, by (19), gives 


^7fcjln(nj/z^j) < \n{vk/vk), 
j 

YQkjt^jMvj/l^j) < l^k^n{vk/l^k), 

j 


In 



QkjVj 

j 


hence £ < £ (as Qkj = !)• Thus, the entropy is a non-decreasing function 
of time in the general case as well. Note that the value of AS" we compute 
here is proportional to Eln(^ where C is the energy gain factor from ([^, i.e., 
we have established our claim ([^. 


An intuitive explanation for the entropy growth at an adiabatic change of 
parameters of a non-ergodic system can be obtained as follows. Let us think 
of an adiabatically changing Hamiltonian system as a gas of non-interacting 
particles (different particles correspond to different initial conditions). As there 
is absolutely no interaction, there is no equilibrium distribution in energies. 
However, in the ergodic case one still recovers the conservation of entropy. In 
the non-ergodic case, the particles are, at each value of the parameter r, in 
different macro-states which correspond to different ergodic measures Hir over 
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which the fast dynamics is averaged. Thus, our gas can be considered as a 
mixture of different phases or fractions; the adiabatic change of parameters 
can lead to particles changing their phase, so the relative densities of each 
fraction in the gas can vary, and this naturally leads to the entropy growth. 


3 Numerical Simulations 

3.1 Quartic Potential 

We verify our model in numerical experiments performed with 

H{p, q, Tit)) = y + y + ^ (gf + + ^qlql, (22) 

where q = ( 91 ,^ 2 )) P = iPi-,P 2 )-, and r(t) = (a(t),6(t)). As the potential 
is homogeneous, the Hamiltonian has scaling properties, i.e., the dynamics 
of the frozen system for any positive energy E can be determined by simple 
scaling from the dynamics on the energy shell E = \. For frozen values of the 
parameters, this system has been thoroughly studied [39lll0]. The system is 
integrable in two situations: 

1. For a = b the conserved quantity as an angular momentum; 

2. For 5 = 0 the two degrees of freedom are uncoupled. 

There is also a chaotic regime; for example, for a = 0.01, 5=1 the system 
exhibits exponential decay of correlations. Thus, we can change a and 5 in such 
a way that the system will undergo a transition between chaotic and integrable 
regimes, see Fig. 


Table 1: Parameter Values used in the simulations. 


Quantity 

Meaning 

Value 

T 

period of oscillation of parameters 

400 

h 

Integration step 

10-^ 

N 

number of points in the ensemble 

N = 2 xW 

Eq 

initial energy 

3 X 10^ 


In the numerical simulations, we always start with initial conditions uni¬ 
formly distributed at the energy level Eq. Throughout the paper (•) stands for 
the ensemble average with respect to these initial conditions. All the relevant 
parameters for the simulations and their physical meaning can be found in 
Table [TJ 
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Figure 3: Parameter plane (a, b) for the quartic Hamiltonian and snapshots 
of dynamical behaviour. The insets show the dynamics of the Poincare map on 
the cross-section y = 0 {at Py > 0). Different choices of the parameters a and b can 
lead to integrable and chaotic dynamics. 

3.2 Symplectic integrator 

We implemented an explicit fourth-order method for the integration of equa¬ 
tions of motion. This method was used to eliminate non-symplectic effects 
while providing the accuracy of a fourth-order integration step. Our Hamilto¬ 
nian is separable, meaning that it can be written in the form 

H{p,q) = A{p) + V{q), 

with A being the kinetic energy and V the potential energy. Introducing the 
symbol 2 = {q,p) and the operator, Dh = {-,11}, which returns a Poisson 
bracket of the operand, the Hamilton’s equation can be written as i = Dhz. 

The formal solution of this set of equations is given as z{h) = exp{hDH)z{0). 

For this class of Hamiltonians the solution reads 

z{h) = exp[h{DA + Dv)]z{0). 

The sympletic integration scheme approximates the time-evolution opera- 
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tor exp[/i(DA + Dy)] in the formal solution by a product of operators as 

k 

ex.p[h{DA + Dv)] = exp{cihDA) ex.p{dihDv) + 

i=l 

For the fourth order integrator m we have 

1 , 1 - 2^/3 

Cl = C4 = — --7:7- and 02 = 03 = —--77— 

2(2 - 2V3) 2(2 - 2V3) 

1 2 V 2 

= ^3 = ^_ 21/2 ’ "^2 = - and di = 0 . 

(23) 

As our Hamiltonian is non-autonomous, we therefore we can introduce pt, 
a canonically conjugate variable of t, so the extended Hamiltonian becomes 

^{P, Q, t) = A{p) +pt + V{q, t). 

This quantity must be conserved by the integration scheme up to the given 
accuracy. We checked the accuracy both with the time-independent and time- 
dependent Hamiltonian and found excellent agreement with the order of the 
method, that is, the Hamiltonian "K conserved up to order 0{h^) per integra¬ 
tion step. Thus, the changes of the energy we observe in the time-dependent 
Hamiltonian system are due to dynamics and are not induced by the numerical 
scheme. 


3.3 Integrable Regime 

If the parameters a and b are periodic functions of time, when keeping the 
Hamiltonian within a given regime of integrability we have observed no energy 
growth. For example, choosing the following parameters 


b = 0.0 and a 


A 


1.1 -h cos ( — 


the two systems are decoupled with A > 1. For example, choosing A = 1 a 
typical behaviour of the energy can be seen in Figure]^ We have experimented 
with A ranging from 1 to 50 we found that no energy growth. Another possibil¬ 
ity is a = 5 depending periodically on time with period T. Within this regime 
we found no energy growth. For example, choosing b = Ax [1.1 -|-cos(27rt/r)]. 
We have experimented with A ranging from 1 to 50 and found no energy 
growth. 
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Figure 4: Lack of energy growth in the integrable regime. We change param¬ 
eters so that the frozen system is integrable for all values of the parameters. In these 
simulations we chose A = 1. The ensemble averaged energy versus time is shown. 


3.4 Strong chaos and polynomial energy growth 

If we change parameters in such a way that for each frozen value of the param¬ 
eters the Hamiltonian remains strongly chaotic, we observe only a slow energy 
growth. For instance, for the parameters 


a = 0.01 and 6(t) = 1.5-|-cos(27rt/T) 


the ensemble energy growth (En/Eo) versus the number of periods n is shown 
in Fig. which clearly reveals a quadratic growth of the ensemble energy. For 
each initial condition, we also measure the energy gain after n periods: 


En 

Eo 


n—1 


n 


Ek 

Ek-i 


and compute the exponential rates 

r(n) = ^ In [En/Eo] . 

In the inset of Fig. we plot the exponential rate r for two distinct initial 
conditions. As we see, r converges to zero corroborating the lack of exponential 
energy growth. 


3.5 Mixed phase Space and Exponential Accelera¬ 
tion 

The most of our numerical experiments correspond to the case where the pa¬ 
rameters visit both chaotic and integrable regions in the parameter space. We 
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Figure 5: Polynomial energy growth in the ergodic regime. We change pa¬ 
rameters so that the frozen system remains chaotic, with no visible stability islands. 
The ensemble averaged energy versus time is shown. The solid (red) line is the least- 
square £t by a quadratic polynomial. In the inset the rates r{n) = ^ In [En/E q] are 
shown for two trajectories for a larger number of periods. Both rates tend to zero, 
corroborating the lack of exponential acceleration. 


demonstrate that in this case the acceleration process is in a good agreement 
with the Geometric Brownian Motion (GBM) model, as predicted by the the¬ 
oretical arguments in Sections o and The parameters are changed along 
the cycle displayed in Fig. which is described by 


a = 




\ T J 
0 - 0 , 


otherwise 


(24) 


together with 



if sin (^) > 0 
otherwise 


(25) 


In Figs. [^and[^ we show results for oq = 0.1, A = 1. 

The exponential energy growth is clearly seen in Fig. Note that two 
distinct growth rates are observed. This is a characteristic signature of mul¬ 
tiplicative random processes: when the averaging is performed over a finite 
ensemble, process Q yields 


the rate p"*" = InE^" for small n, then a crossover to the lower rate p = Eln(( for large n 


(typically, the rate p is established only after very large n; e.g. it was unachiev¬ 
able in our experiments, so we can see only the beginning of the cross-over pro¬ 
cess) . The cross-over phenomenon is well known [38t l3T| [22] ; it is explained by 
the fact that the variance Var(£'„) grows much faster than E(£'„). In partic¬ 
ular, when process ([^ is GBM (i.e., InC is normally distributed) the variance 
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Figure 6: Exponential energy growth. We show the ensemble energy in log- 
scale versus the number of parameter cycles. The parameters cycle described by 
equations and for oq = 0.1, A = 1, see Figure]^ For these parameters, the 
frozen Hamiltonian exhibits chaotic dynamics (along a = 0.01), quasiperiodic motion 
(along 6 = 0), and mixed behaviour (along the connecting arc). Repeating the cycle 
multiple times, we observe exponential growth of energy. For the hrst 70 cycles, 
we observe exponential growth at a higher ensemble rate p+, then a crossover to a 
lower rate starts. This crossover is a characteristic feature of a Geometric Brownian 
Motion. 


is given by 

Var(E„) = - l) (26) 

where 

= 2(p+ - p) (27) 

In Fig. [^we show the energy distribution after few periods of parameter 
oscillations (initial distribution is a Dirac delta-function, i.e., we start with 
uniformly distributed initial condition on the energy shell Eq). One can see 
that the logarithm of the energy hiEn is well described by a Gaussian distri¬ 
bution, already after the first cycle. This gives a qualitative agreement with 
our claim 0. 

In order to make a quantitative check of our GBM model, we estimate the 
parameters of the GBM from the data. The values of ensemble rate p+, single 
particle rate p and the variance are estimated as 

/5=-/ln^y p+ = -ln/^y and On = \sx{En/ Eq). 
n\ Eq / n \Eq / 
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Figure 7: Distribution of the logarithm of the energy. Initially all particles 
start at the same energy level Eq, then after n cycles of the parameters, that is, after 
a time riT, we measure the logarithmic energy gain \\i{En/E q). The distribution of 
the energy gains v is shown at n = 1 in inset a), n = 3 inset b), and n = 5 for inset 
c). The full line (red) shows the Gaussian fit. 


In our experiment, the values of the estimators stabilise already at the first 
cycle. Next, we estimate a relative error with which the GBM predictions 
p^,([27l) hold: 

q2 _ g2p+n J'g2(p+-p)n _ 


X = 


02 


In Table 2, we show the estimators quantities for the first five cycles. The 
estimators give excellent agreement with the GBM model, with the relative 
error y about 3%. 


Table 2: Value of the estimators for the GBM for parameters in (24) and (25) for 
Oo = 0.01 and A = 1. 


n 


P 

02 

X 

1 

0.0391 

0.0256 

0.0305 

0.029 

3 

0.039 

0.0234 

0.1175 

0.055 

5 

0.0387 

0.0238 

0.2307 

0.025 


In order to see how sensitive is the GBM prediction to the variations of 
parameters of the system, we use the following procedure. From the data, 
using the values obtained for 0, we estimate cr^ by formula (26). Denote this 
estimated value Then we can systematically check the GBM relation (27), 
which is a consequence of 0. To this end, we notice that by changing the 
parameter oq we can tune the ensemble rate p"*". In FigJ^ we vary oq in 
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the interval 10“^] and compute the ensemble rate after one cycle of 

parameters oscillation. As oq decreases, the volume of the phase space (of the 
frozen system) increases, enhancing the exponential rate p'^. Then, for each 
value of oo we compute p'^, p and a and verify relation (27). The result can 
be seen in Fig. [8>), which shows excellent agreement with the GBM model. 




Figure 8: Geometric Brownian motion for various parameter values. Inset 
a) shows the estimated ensemble rate versus the parameter oq. For each value of 
Oo we estimate the GBM parameters and check relation 27 In inset b) we plot the 
estimate p obtained from data against the parameters p^ and also obtained from 
the data. The full circle shows the results whereas the full line (red) is the identity 
line. The results are in excellent agreement with the GBM model. 
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